function [vp, vr, vu] = fn_iteration(n,h0,vk,vFT,ck,cX,Num_point)

a = fn_inveta(n,h0);
f = @(x) fn_f(n, x, a, vk, vFT, ck, cX); 

p_left = 0;
p_right = 1000;

% vp = [p_left, p_right];
% [r_left, u_left] = f(p_left);
% [r_right, u_right] = f(p_right);
% vr = [r_left, r_right];
% vu = [u_left, u_right];

vp = linspace(p_left, p_right, 100);
vr = zeros(1, length(vp));
vu = zeros(1, length(vp));

for i = 1:length(vp)
    p = vp(i);
    [r, u] = f(p);
    vr(i) = r;
    vu(i) = u;
end

while 1
    
    vd = abs(vr(2:end) - vr(1:end-1));
    
    [~, max_region] = max(vd);
    
    ind_left = max_region;
    ind_right = max_region + 1;
    
    p_left = vp(ind_left);
    p_right = vp(ind_right);
    p_mid = 1/2*(p_left + p_right);
    
    [new_r, new_u] = f(p_mid);
    
    vp = [vp(1:ind_left), p_mid, vp(ind_right:end)];
    vr = [vr(1:ind_left), new_r, vr(ind_right:end)];
    vu = [vu(1:ind_left), new_u, vu(ind_right:end)];

    if length(vp) >= Num_point
        break;
    end

end

end